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A prolonged timing of millisecond pulsars has revealed low-frequency uncorrelated (in¬ 
frared) noise, presumably of astrophysical origin, in the pulse arrival time (PAT) residuals 
for some of them. Currently available pulsar timing methods allow the statistical parameters 
of this noise to be reliably measured by decomposing the PAT residual function into orthog¬ 
onal Fourier harmonics. In most cases, pulsars in globular clusters show a low-frequency 
modulation of their rotational phase and spin rate. The relativistic time delay of the pulsar 
signal in the curved space time of randomly distributed and moving globular cluster stars 
(the Shapiro effect) is suggested as a possible cause of this modulation. Extremely important 
(from an astrophysical point of view) information about the structure of the globular cluster 
core, which is inaccessible to study by other observational methods, could be obtained by 
analyzing the spectral parameters of the low-frequency noise caused by the Shapiro effect 
and attributable to the random passages of stars near the line of sight to the pulsar. Given 
the smallness of the aberration corrections that arise from the nonstationarity of the gravi¬ 
tational field of the randomly distributed ensemble of stars under consideration, a formula 
is derived for the Shapiro effect for a pulsar in a globular cluster. The derived formula is 
used to calculate the autocorrelation function of the low-frequency pulsar noise, the slope 
of its power spectrum, and the behavior of the cr z statistic that characterizes the spectral 
properties of this noise in the form of a time function. The Shapiro effect under discussion 
is shown to manifest itself for large impact parameters as a low-frequency noise of the pul¬ 
sar spin rate with a spectral index of n = —1.8 that depends weakly on the specific model 
distribution of stars in the globular cluster. For small impact parameters, the spectral index 
of the noise is n = —1.5. © 2006 Pleiades Publishing Inc. 
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INTRODUCTION 


Explaining the astrophysical origin of the low frequency uncorrelated pulsar timing noise 
is an extremely important, but challenging problem of modern pulsar astronomy. It is 
well known that millisecond and binary pulsars can be used as the most stable natural 
frequency standards (Davis et al. 1985; Rawley et ah 1987; Kaspi et al. 1994; Petit 
and Tavella 1996; Ilyasov et al. 1998; Kopcikin 1999). However, obtaining an unbiased 
and statistically significant estimate of their stability depends considerably on the accurate 
determination of the parameters of the random low-frequency noise that is detected in the 
pulse arrival time (PAT) residuals observed on long time scales. In most cases, the low- 
frequency (infrared) noise can be described by an additive set of components each of which 
is described by a simple model of the autocovariance function with a power-law spectrum, 
S(f) = A n f~ n , where the index n takes on integer values from unity and higher (Kopcikin 
1997). Noise of different physical natures corresponds to different spectral indices. For this 
reason, studying the spectral parameters of the infrared noise is an important problem of 
observational astrophysics. It enriches significantly the metrological studies of the pulsar 
time scale and encompasses many divisions of modern astrophysics, the most important of 
which are cosmology and gravitational wave astronomy (Kopcikin 1997a; Lorimcr 2001). 

A self-consistent approach to estimating the stability of the rotational and orbital param¬ 
eters of pulsars and pulsar noise models has been developed by many authors (Groth 1975; 
Cordes 1980; Bertotti et al. 1983; Blandford et al. 1984; Taylor 1991; Kopcikin 1997a, 1997b, 
1999; Ilyasov et al. 1998; Kopcikin and Potapov 2004). They solved the main theoretical 
questions of the analytical observational data processing technique in the time domain. The 
data processing in the time domain is more informative than that in the frequency domain, 
since both the stationary and nonstationary noise components are taken into account in the 
former case, while an adequate data analysis in the frequency domain and obtaining unbiased 
estimates of the measured parameters of the time series are possible only for a stationary 
random process. 

It should be noted that the accuracy of determining the pulsar PATs is currently very high 
and approaches ~ 50 ns on a time scale of several years (Bailes 2003). Therefore, we need a 
high-precision model for processing the PATs that are affected by many factors, such as the 
Earth’s orbital and rotational motions, the proper motion of the pulsar, the gravitational 
potentials of the Solar system and other gravitating objects at the point of observation and 
along the pulse propagation path, and the spatial distribution of the interplanetary and 
interstellar media. The classical and relativistic effects related to the orbital motion of the 
pulsar around the barycenter of the binary and the proper motion of the latter in space, the 
propagation of the radio emission from the pulsar in the gravitational held and atmosphere 
of its companion, and the precession of the pulsar spin axis, which causes the pulsar pulse 
shape to change with time, are added for pulsars in binary systems. 

Based on the relativistic theory of astronomical reference frames and time scales (Kopcikin 
1988, 1989a, 1989b; Brumberg and Kopeikin 1989, 1990), whose improved and extended ver¬ 
sion was taken by the General Assembly of IAU-2000 as the basis for relativistic astronomical 
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algorithms and was described in detail in review papers (Soffel et al. 2003; Kopcikin and 
Vlasov 2004), Doroshenko and Kopcikin (1990, 1995) developed a pulsar timing algorithm 
that includes the above effects, to within 10 ns, and that is suitable for processing the obser¬ 
vations of both single and binary pulsars. This theory was effectively used to systematically 
develop the TIMAPR software package (Doroshenko 1997), which is designed for pulsar data 
processing (Larchenkova and Doroshenko 1995). The independent TEMPO code was devel¬ 
oped at the Princeton University, USA (Taylor and Weisberg 1989). Both codes are widely 
used by Russian and foreign scientists at various radio observatories worldwide (Lorimer 
2001). 

The standard procedure for estimating the pulsar parameters on fairly short time scales 
is based on the assumption that white noise dominates in the PAT residuals. However, a 
fairly long monitoring reveals components of the correlated infrared noise of astrophysical 
origin almost in all pulsars, whose spectrum differs from the white noise spectrum (Cordes 
and Downs 1985; D’Alessandro et al. 2001). The spectrum of this correlated (infrared) noise 
diverges at zero frequency. This infrared catastrophe compels the researches to reconsider 
the standard methods of spectral analysis and forces them to resort to various kinds of 
regularization procedures that allow the divergence to be avoided when modeling the noise 
spectrum (Kopeikin and Potapov 2004). 

One of the most important problems in modern pulsar timing is to separate out the 
infrared noise in the PAT residual spectra and to determine its amplitude and spectral index. 
Although this is a challenging problem, its solution can provide substantial information 
about the physical processes inside neutron stars and in the interstellar medium in the 
path of pulsar pulse propagation and help to detect low-frequency gravitational waves and 
other, no less interesting gravitational effects. The most suitable objects for solving this 
problem are millisecond pulsars with a very high spin stability and, as a result, a low level 
of intrinsic rotational noise (Guinot and Petit 1991). Millisecond pulsars were discovered in 
various regions of our Galaxy. However, the population of millisecond pulsars in globular 
star clusters is currently most representative. The first millisecond pulsar was discovered in 
1987 in the core of the globular cluster M28 (Lyne et al. 1987). A systematic search for 
pulsars both in M28 and in other globular star clusters began after this discovery. Twenty 
four globular clusters in which more than 80 pulsars have been discovered are known to date 
(Freire 2004). The globular clusters 47 Tuc (NGC 104) and M15 (NGC 7078), which contain, 
respectively, 22 and 8 pulsars (Taylor et al. 1993; Camilo et al. 2000; Freire 2004), and the 
globular cluster Terzan 5, in which 26 pulsars were discovered (Ransom et al. 2005), are the 
record-holders in the number of discovered pulsars. 

If the inherent causes of the pulsar spin instability are ignored, then three physical causes 
of the low frequency timing noise produced by external effects can be suggested. First, the 
stochastic Shapiro effect, i.e., the random spread in the total time it takes for the pulsar radio 
signal to pass through the fluctuating gravitational held of the globular cluster stars that 
are distributed and move randomly in a certain vicinity of the line of sight to the pulsar; 
second, the stochastic gravitational perturbations in the pulsar velocity and acceleration 
produced by close passages of globular cluster stars near the pulsar itself (Joshi and Rasio 



1997; Rodin 2000); and, third, the random fluctuations of the interstellar medium (Smirnova 
and Shishov 2001). In this paper, we concentrate on detailed calculations and analysis of 
the spectrum of the pulsar noise produced by the stochastic Shapiro effect. The amplitude 
of this noise depends significantly on the radial star density distribution inside the cluster, 
which is unknown and model-dependent in most cases. From general considerations, we may 
expect the noise amplitude to be comparable to the gravitational radius of a typical cluster 
star, i.e., ~ 10 /xs. However, this value may be an order of magnitude higher, since the 
noise from the Shapiro effect has a cumulative property and is directly proportional to the 
number of stars on the line of sight to the pulsar that produce this noise. There is no doubt 
that the increase in the noise amplitude correlates with the characteristic time scale of its 
manifestation. 

The random process produced by star passages should be considered as a special case of 
gravitational lensing. When applied to microarcsecond astrometry, this case was considered 
by Sazhin et al. (1998, 2001). In this paper, we calculate the noise spectrum, as applied to 
pulsar observations. Knowledge of the theoretical power spectrum of the stochastic Shapiro 
effect for a pulsar in a globular cluster will allow us, in the case of its direct measurement, to 
obtain very important astrophysical information about the structure of the globular cluster 
core, which is inaccessible to other observational methods, and to analytically extend the 
Salpeter mass function for globular clusters toward the low-mass stars comparable in mass 
to Jupiter. Observation of the stochastic Shapiro effect will also allow the mass of the dark 
matter that is possibly concentrated near the globular cluster cores to be estimated. 

Kopeikin and Schafer (1999) were the first to derive an exact formula for the Shapiro effect 
produced by a moving massive body and showed that the positions of the light-deflecting 
gravitating bodies (stars) should be taken not at the pulsar pulse arrival time to the ob¬ 
server, but at the corresponding delayed time due to the finite speed of propagation of the 
gravitational perturbation (Kopeikin 2001). This relativistic effect has received convincing 
experimental confirmation through high-precision VLBI measurements of the deflection of 
light from a quasar by the gravitational field of a moving Jupiter (Famalont and Kopeikin 
2003). We took into account this effect when deriving an expression for the total relativistic 
delay time of the pulsar signal in the gravitational field of arbitrarily moving globular cluster 
stars. 

This paper is structured as follows. In the next section, we briefly consider the pulsar tim¬ 
ing model and the procedure for estimating its parameters in the presence of low-frequency 
noise and introduce the concept of a z statistic, which serves as a quantitative measure of 
pulsar instability. Subsequently, we write out a formula for the cumulative Shapiro effect 
for a pulsar in a globular cluster and an expression for the cr z statistic in the case where 
the physically significant noise is attributable to the Shapiro effect from random passages 
of stars near the line of sight. In conclusion, we use the a z statistic to estimate the noise 
spectrum and discuss prospects for numerically analyzing this effect. 
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PULSAR TIMING, ESTIMATION OF PULSAR PARAMETERS, AND THE 

g z STATISTIC 


Let the observations begin at time t 0 - The rotational phase of a pulsar is specified by a 
time polynomial: 

iV(f) = is p ^s + —UpQ 2 + —is p Q 3 + — LpiJ 4 + UpfoC^s) + 0(A 5 ) (1) 

where S = ^s(t) is the pulse emission time in the proper time scale of the pulsar; t is the 
pulse arrival time in the barycentric time scale of the Solar system; is p , u p , is p and 'is p are the 
pulsar spin rate and its derivatives taken at time S = 0; and 02 (9) is the intrinsic noise of 
the pulsar rotational phase, spin rate, and its time derivatives. 

Note that the pulsar proper time 9 is not directly observable. The barycentric time, t, 
which is related to the proper time of the observer (atomic time) by a relativistic trans¬ 
formation (Kopeikin 1989b; Brumberg and Kopeikin 1990; Kopeikin and Vlasov 2004), is a 
measurable quantity. The relationship between the pulsar proper time S and the barycentric 
time t is established by solving the equations for the light geodesics that describe the prop¬ 
agation of the radio pulse from the pulsar to the observer (Kopeikin 1990; Doroshenko and 
Kopeikin 1990). This relationship allows the observed pulsar rotational phase N(t) to be 
expressed as a function of the barycentric time t (Doroshenko and Kopeikin 1990; Kopeikin 
1999): 


1 1 

N (t) = Nq + ist + —ist 2 + —ist 3 + is£(t), (2) 

2 o 

where N 0 is the initial rotational phase of the pulsar; is, is and is are the barycentric 
pulsar spin rate and its derivatives at the initial epoch of observations T; and e{t) is the 
total, physically significant additive noise of the pulsar rotational phase, with the intrinsic 
noise of the pulsar 02(f) being one of its components. We emphasize that the noise e(t) is 
produced by both external and internal factors. By the external factors we mean all those 
factors that affect the propagation of the pulsar radio pulse in a random way. These also 
include the noise from the passages of globular cluster stars, which will be considered in 
detail in the next sections. The internal factors responsible for the pulsar noise are related 
to the pulsar spin mechanism and are not considered here. Lyne and Graham-Smith (2004) 
gave a comprehensive overview of the causes of the pulsar spin instability. 

The PAT residual r(t) is the difference between the observed pulsar rotational phase, N obs , 
and its theoretical value, N(t,9), predicted by the timing model and specified by Eq. (2) 
divided by the pulsar rotation rate, is: 

IS 

where 9 = {9 a ,a = 1, 2, k} denotes a set k of measurable parameters {k = 4 in the 
model represented by Eq. (2)). 
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Let us introduce the quantities (5 a = 0* a — 9 a that are the corrections to the unknown true 
values of the parameters 6 a and the timing model fits 


^a(t) 


d N 

df n 


e=e* 




where 8* are the least-squares estimates of the parameters (Kopeikin 1999). The parameters 
and fits are given in the table. 

Assuming that m equally spaced and comparable (in accuracy) pulse arrival times are 
measured for N complete rotations of the pulsar around its axis, we have mN residuals: 
r\ = r(ti), where i = l,2,...,miV. The least-squares method yields the parameters /3 a (r) 
(Bard 1974): 

4 mN 

Pair) = L abMti)£(ti), a = 1, 2 , 3,4, 

6=1 i= 1 


where the information matrix 

mN 

L a b(r) = 

i =1 

r is the total observing time. The correlation matrix of the parameters is defined as M ah = 
{PaPb), where the angular brackets denote an ensemble average. 

The behavior of the PAT residuals can be described best using the so-called cr~ statistic. 
It is defined as the weighted root-mean-square value of the coefficients of the cubic terms 
in the time polynomial fitted to the observed pulsar phase divided into the segments that 
correspond to equal time intervals r (Matsakis et al. 1997). The cr 2 statistic is formally 
defined (Matsakis et al. 1997) as 

a ’ {T) = 

where M 44 is a diagonal element of the correlation matrix M a b. This diagonal element is 
defined by 


c= 1 d= 1 


mN mN 




L*=i 3 = 1 


( 3 ) 


where ) = (e(L)e(tj)) is the autocorrelation function of the random process e{t). 

For a stationary noise process, the autocorrelation function depends not on the individual 
times ti and tj , but only on their difference r: 

r) = (e(t + r)e(t)) = 3?(r), 


where t — tj, r — U — tj. 
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List of basic functions and parameters used to fit the parameters in the pulsar timing 

model specified by Eq.(2) 


Parameter 

Fit 

(5i = 5N 0 /v 

V'i (t) = i 

(3 2 = 5v/v 

^2 (t) = t 

[3 3 = bvfv 

^3 (t) = t 2 

04 = 5v/v 

V>4 (t) = t 3 


THE SHAPIRO EFFECT FOR A PULSAR IN A GLOBULAR CLUSTER 

The relativistic time delay of an electromagnetic signal in a static, spherically symmetric 
gravitational held of a point mass is called the Shapiro effect (Shapiro 1964). A gener¬ 
alization of the formula for the Shapiro effect for the propagation of light in a variable 
gravitational held of an arbitrarily moving body was found by Kopeikin and Schafer (1999) 
and experimentally confirmed by Fomalont and Kopeikin (2003). In this paper, we use the 
Kopeikin-Schafer formula. 


Formula for the Generalized Shapiro Effect 


The expression for the travel time of an electromagnetic signal between fixed points, 
(t o,Xq) and (t,x), in the gravitational held of an arbitrarily moving body (star) derived by 
Kopeikin and Schafer (1999) is 

t ~ t 0 = \x- x 0 \ + A(Mo), 

where \x — x$\ is the coordinate distance in the background Euclidean space between the two 
points, xo and x; and A(f,f 0 ) is the relativistic time delay attributable to the gravitational 
held of the moving bodies. 


For a pulsar in a globular cluster, the Shapiro delay produced by the gravitational held 
of the moving cluster stars can be written as (Kopeikin and Schafer 1999) 


A(Mq) 


N 


2GM a ^ r a 
* C 3 r 

a= 1 r 0 a 


(k 0 f a ) 

( k 0 f 0a ) ’ 


(4) 


f a = X(t) ~ X a (s), 

f 0a = X 0 (t 0 ) - f a (S 0 ), 

where we discarded the terms that are proportional to the velocity of the bodies and that 
appear in the amplitude of the logarithmic function, because they are small. The coordinates 
of the gravitating bodies are calculated at the delayed times s and so, which are defined by 
the equations for isotropic gravitational held characteristics: 


t = s+ | x(t) - x a (s) |, 
h) — s o + \x 0 (t 0 ) — f o (s 0 )|. 


( 5 ) 
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Here, to is the photon emission time, x specifies the observer’s position relative to the 
globular cluster barycenter (GCB) at the observation time of the pulsar radio pulse, x a 
specifies the position of star a relative to the GCB, xo specifies the pulsar position relative to 
the GCB at the radio pulse emission time, and the unit vector ko specifies the direction of the 
rectilinear, gravitationally unperturbed radio pulse motion from the pulsar to the observer 
and is defined by 



x-x 0 (t 0 ) 

x-x 0 (t 0 )\' 


Let us expand r a , r a , and ko in a Taylor power series of x a /R -C 1, where R = \x(t)\ is 
the distance from the observer to the GCB, much larger than the cluster size: 

k 0 — K — K x (f x K) + 0(g), 

Ta = R- R(Ma) + \R(L X Kf + 0(g), 
where K = x{t)/\x{t)\ is a unit vector, f = ^§1* and g = g|j. 

Let T be a fixed start time of observations. Let us expand x(t) and x a (s) in a Taylor power 
series of (t—T), provided that the following conditions are satisfied: \x(T)\ S> \v(T)(t—T)\ S> 
\v(T) (t—T) 2 1. Retaining only the linear (in velocities) terms and using the following notation: 
Ro = \x(T)\, Ko = Ro/Ro is a unit vector, v is the observer’s velocity, v a is the velocity of 
star a, and vq is the pulsar velocity, we then obtain 


K = IK + 


x(t) = x(T) + v(T)(t-T) + ..., 
x a (s) = x a (S) + v a (s - S) + ..., 
R = R 0 + (K 0 v(T))(t-T) + ..., 
v(T) - K 0 (K 0 v(T ))' 


Ro 


(t — T) + 


where S is the delayed time related to the start time of observations T by the gravitational 
cone equations (5). 

Using the equations for gravitational field characteristics (5), we obtain the following 
expression for the numerator of the fraction in the argument of the logarithm in Eq. (4) 
after algebraic transformations: 



x , o d a v a (S) 

da%a(S^) 


(t-T) + ..., 


where d a = (K 0 x j^) x K 0 is the impact parameter of the radio pulse trajectory with respect 
to the GCB. 


Acting in a similar way, we expand fo a and ro a in a Taylor series, retaining only the 
linear (in velocity) terms, and obtain an expression for the denominator of the fraction 
in the argument of the logarithm in Eq. (4) after simple transformations. Finally, once 
all terms have been reduced to the same time To, which defines the pulsar pulse emission 



time corresponding to the start of observations T, we obtain the ultimate expression for the 
random process produced by the relativistic delay during the propagation of a radio pulse 
in a variable gravitational held of moving globular cluster stars: 


2GM a 

a= 1 


C 


In < 1 + 2 (t — T) 


dq(T 0 )v a (T 0 ) 
d a (T 0 )x a (T 0 ) 


In < 1 + 


e{t) = A(f, t 0 ) — A(T, T 0 ) 
-q~ K 0 )V 0a (T 0 ) 


Qoa + dioQo a 


<t-T) 


( 6 ) 


where V 0a = v 0 (T 0 ) - t? a (T 0 ), Q 0 a = AiC^o) - T 0 (T 0 ), Q 0a = (Q 0a Qoa) 1/2 , Q = Qoa/Qoa is a 
unit vector. 


Autocorrelation Function of the Shapiro Effect 


Let us pass to a new model of the random process e{t) defined by the formula 

e{t) = e(t) - (e(t)), 

where the angular brackets denote a statistical ensemble average, and we assume that the 
mean value is (e(t)) = 0. The autocorrelation function in Eq. (3) can then be written as 

t) = j dm a dx a dv a f(m a , x a , v a )e(t, m a , x a , v a )e(t + r, m a , x a , v a ), (7) 

where we assume that the statistical ensemble of stars is defined by uncorrelated parameters, 
so the distribution function can be fitted by the product of three statistically independent 
distribution functions: 

f(m a ,x a ,v a ) = Af(m a )f(x a )f(v a ); 
we find the normalization numerical coefficient A from the condition 

A = fdm^U = 


No integration limits are specihed in Eq. (7), but we assume that they are known (see below) 
and specify the range of statistical ensemble parameters. 

A globular cluster is characterized by two quantities: the cluster core radius r c , which is 
defined as the distance at which the surface brightness is half its central value, and the tidal 
radius r t , at which the surface brightness is zero. Given the mass of the Galaxy Mq , the 
cluster mass M c , and the Galactocentric distance of the cluster center Rq, the tidal radius 
can be calculated as follows: rj = j^Rg- Equation (7) is integrated in the following limits: 
m a = [O.IMq — 1OM 0 ], | va = [0 — 4cr] and \x a \ = [0 — r t ], where a is the stellar velocity 
dispersion of the cluster, and M 0 is the solar mass. 

We use the Salpeter function f(m a ) ~ m~ 2 35 as the mass function of globular cluster 
stars. Let us consider two model density distributions of a globular cluster: the model of 
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an isothermal sphere with a core and the King model. In spherical coordinates, the density 
distribution for the model of an isothermal sphere with a core is (Spitzer 1990) 


f( r a) 


Po r l 

r 2 i r 2 ’ 

' c ' a 


( 8 ) 


where p is the core density of the globular cluster. The velocity distribution in the model of 
an isothermal sphere has a Gaussian (Maxwellian) profile, 


f(Va) 


1 

2rr a 2 


exp 


+ v 2 ay + V 2 az 
2a 2 


(9) 


It is the Maxwellian distribution that is typical of unrelaxed systems, which the globular 
clusters in our Galaxy are. However, the model of an isothermal sphere is not devoid of 
shortcomings. The main shortcoming is the assumption of its infinite radius in geomet¬ 
ric space and velocity space, which is not physically justified. This shortcoming can be 
eliminated by ”truncating” the argument of the velocity distribution function (9) by the 
characteristic star escape velocity v e . This allows the isothermal sphere to be made finite in 
velocity space while preserving the isotropic Maxwellian distribution. From physical consid¬ 
erations in the velocity range v a > v e , the distribution function f(v a ) must be close to zero 
and can be fitted by a truncated Maxwellian function: 

e -v 2 J ° 2 _ e -vl/a ‘ 2 

1 _ e -w 2 G 2 ’ 


f(Va) = 


where the normalization was chosen in such a way that f(v a ) = 1 at v a = 0. The models of 
star clusters based on the truncated Maxwellian distribution were calculated by King (1966) 
and are called the King models. The density distribution in the King models is also bounded 
in space by the tidal radius r t and is defined as (King 1966) 


(1 + T 2 ) arccos \ f- - a/T 2 - (r a /r c ) 2 y/l + (r a /r c ) 2 

f(r a ) = Po -- X - i== -n-, (10) 


[1 + (r a /r c ) 2 ] 3/2 (1 + T 2 ) arccos J ^2 - T 


where T = r t /r c . For the subsequent calculations, we use the globular cluster 47 Tucanae, 
in which the largest number of pulsars have been discovered to date. This cluster has the 
following parameters: r c = 0.52 pc, p 0 = 6 x 10 4 M o /pc 3 , r t = 60.3 pc, and a = 10 km s _1 ; 
the distance to the cluster center is Rq = 4.1 kpc (Harris 1996). To determine the slope of 
the power spectrum for the noise process under study and the behavior of the a z statistic, 
we must calculate the integral in Eq. (7). Analysis indicates that this integral cannot 
be calculated analytically in general form and the effect under study must be simulated 
numerically, which is done here. The random process attributable to the Shapiro effect can 
be estimated analytically if we restrict our analysis to star passages far from the line of 
sight to the pulsar. The random process under study can be separated into long- and short- 
period parts. The long-period part of the process corresponds to the case where the impact 
parameter da of the pulsar pulse trajectory is large for the cluster stars (analytical case). 
The short-period part of the process corresponds to the case of small impact parameters. 
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Large Impact Parameters 

Let us make the simplifying assumption that the distributions of globular cluster stars 
in space, velocity, and mass are uncorrelated. This simplifies significantly the form of the 
distribution function in Eq. (7): 

f(rn a ,x a ,Va, rn b ,x b ,v b ) = 5(m a - m b )S(x a - x b )S(v a - v b )f(m a ,x a ,v a ), (a ^ b ) (11) 

where 8(y) is the Dirac delta function. The assumption about large impact parameters of 
the pulsar pulse trajectory relative to the globular cluster stars allows the expression for the 
random process under study ( 6 ) to be expanded in a Taylor power series of (t — T): 

e{t) = a(t - T) + P(t - Tf + 7 (t - T ) 3 + ..., (12) 

where T is the initial epoch of observations and the coefficients a, (3 and 7 are defined by 

2 GM a 2 d a v a (—f/— Ko)Voa 

a = > -^— —-, 

0= 1 C d a X a Qoa + KoQoa 

/j = y 2GM a l x f U-q-Ko)Voa \ " _ / 2d a v a \ 
hi c3 2 \ . Qoa + KoQoa _ V daXa ) 

^ y^ 2GM a l I ( 2d a v a \ 3 U-q-K 0 )V 0 a 

a= l c3 ^ | y d a X a J Qoa + KoQoa 

In what follows, it will suffice to restrict the analysis to only the first term in expansion 
(12), which is linear in time (t — T). Let us write the autocorrelation function for it using 
the expression for the distribution function ( 11 ): 

3fJ(ti,t 2 ) = j dmdxdva 2 (m,x,v) x f(m,x,v)(ti — T)(t 2 — T) — (13) 

J dmdxdva 2 (m, x, v)f(m , x, v) x [tit 2 — T(ti + t 2 ) + T 2 }. 

Let us designate r = \t 2 — ti| and t + = (ti + t 2 ). These designations allow the autocorre¬ 
lation function (13) to be represented as 

S( tl . t2 ) = X (7 - T j-2UT + T>) . (14) 

In this formula, the terms proportional to t 2 + , t + T , T 2 constitute the nonstationary part 
of the noise. The terms proportional to t + contribute only to the initial rotational phase 
of the pulsar, while t+, t + T, T 2 are equal to the products of the fits. According to the 
theorem that was proved by Kopeikin (1999), these products of the fits in the structure 
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of the autocorrelation function do not contribute to the PAT residuals, which depend only 
on the stationary part of the noise. Thus, only the stationary component ^ in Eq. (14) 
contributes to the PAT residuals. For this reason, the random process attributable to the 
regular part of the Shapiro effect for pulsars in globular clusters will manifest itself mainly 
as the noise of the pulsar spin rate. 


Numerical Model for the Pulsar Noise Produced by the Shapiro Effect 

As we noted above, knowledge of the autocorrelation function (7) allows us to determine 
the slope of the power spectrum for the noise process and the behavior of the cr z statistic. 
We calculated the triple integral (7) numerically. We constructed two numerical models for 
the density distribution of the globular cluster 47 Tuc described by Eqs. (8) and (10).We 
also considered two extreme cases: small and large impact parameters. 

Large impact parameters. We calculated the autocorrelation function in two ways: 
(1) in a spherical coordinate system with the origin at the GCB and (2) in a Cartesian 
coordinate system. We equated the volumes of the integration space in both coordinate 
systems and took into account the analytic singularities of the integrand by choosing the 
appropriate path of integration. Indeed, it follows from the form of the expression for the 
Shapiro delay (4) attributable to the gravitational held of the moving cluster stars that a star 
on the line of sight produces an in.nite signal delay. For a numerical model, this necessitates 
cutting out a cylinder of minimum radius r c i in the direction of the unit vector k (i . Given 
the radius of the cutout cylinder, we can easily estimate the cluster mass, which is excluded 
from the analysis. For example, for r r j 1AU, it is ~ 10 4 M q . 

The relative accuracy of calculating the integral was determined by the standard method 
and is 0.9% for the King model and 1% an isothermal sphere with a core. The autocorrelation 
function for the King model and the model of an isothermal sphere with a pulsar at the cluster 
center is plotted in Figs. 1 and 2, respectively. 

Analysis of the plots of the autocorrelation function indicates that it is well fitted by a 
quadratic polynomial for both the King model and the model of an isothermal sphere. For 
example, 5R = 7 x 10 -23 + 1.16 x lCT 20 t + 7 x lCT 23 t 2 for the King model (r = 3 yr) and 
3? = 5 x 10~ 24 + 6.12 x 10 22 7 + 4 x 10 24 t 2 for the model of an isothermal sphere with a core 
(r = 3 yr). As would be expected, our numerical result matches the analytic prediction made 
in the previous subsection and confirms the validity of the numerical integration model. 

A stochastic process can be described both in the time domain by a time-dependent 
quantity h{t) and in the frequency domain by the amplitude of the process H(f) as a function 
of the Fourier frequency /. These two quantities are related by the Fourier transform: 

/ OO 

H(f)e 2mft dt 

O 



( 15 ) 
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According to the Parseval theorem, the total power of the signal will be the same when 
calculated in both the time and frequency domains: 

/ OO /»oo 

\h(t)\ 2 dt= / \H(f)\ 2 df. 

If the function h(t) is real, then we can dehne the onesided spectral power density as 

P{f) = 2\H{f)\ 2 . 

The total power will then be calculated as an integral of P(f) in the Fourier frequency 
limits from 0 to oo. 

To analyze the noise process, we numerically determined the spectral power density P(f) 
and the spectral slope n using the fast Fourier transform, whose algorithm has been well 
developed (Elliott and Ran 1982). Our calculations show that the slopes of the power 
spectrum are n = —(1.78 ± 0.04) for the King model and n = —(1.76 ± 0.04) for the model 
of an isothermal sphere, which are equal within the computational error limits. In Fig. 3, 
logP is plotted against log f for the King model. 

Small impact parameters. The random process attributable to star passages with 
small impact parameters should be considered as a special case of gravitational lcnsing. In 
our case, the characteristic impact parameter is equal to the Einstein-Chwolson radius. For 
a pulsar at the center of a globular star cluster, the Einstein-Chwolson radius is de.ned by 
(Zakharov and Sazhin 1998) 

where D ds is the distance from the pulsar to the gravitating body (star), D d is the distance 
from the observer to the gravitating body, and D s is the distance from the observer to the 
pulsar. For the globular cluster 47 Tuc and a gravitating mass of the order of the solar mass 
M 0 , the characteristic Einstein-Chwolson radius is ~ 1 AU. Therefore, in the direction of 
the pulsar specified by the unit vector ko, we cut out a cylindrical volume with a polar radius 
equal to the Einstein-Chwolson radius. We calculated integral (7) in a cylindrical coordinate 
system with a coordinate grid whose density increased with decreasing polar radius, since the 
stars with the smallest impact parameters make a larger contribution to the noise process. 

The autocorrelation function of the noise process produced by the Shapiro effect is plotted 
in Fig. 4 for impact parameters smaller than 1 AU in the model of an isothermal sphere for 
three observing times r: 1, 3, and 5 yr. The derived time dependences of the autocorrelation 
function cannot be fitted by a quadratic time polynomial, as in the case of large impact 
parameters, because the gravitational held acts on the radio pulse produced by close star 
passages near the line of sight for a short time. The autocorrelation function in Fig. 4 is 
nearly logarithmic and, in all probability, can be interpreted as a dicker noise (Kopcikin 
1997b, 1999). The spectral index of the pulsar noise for small impact parameters is n = 
— (1.55±0.03). Figure 5 shows the slope of the power spectrum for the model of an isothermal 
sphere for impact parameters smaller than 1 AU and an observing time of r = 5 yr. 
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Dependence of the noise on cluster core radius. In conclusion, let us analyze the 
dependence of the behavior of the autocorrelation function for the stochastic Shapiro effect 
on the globular cluster core radius. The numerical model described in the previous section 
was constructed for the globular cluster 47 Tuc whose core radius is r c = 0.52 pc. We 
modeled the behavior of the autocorrelation function for three different cluster core radii, 
r c = 0.08, 0.1, and 0.52 pc. The results of our calculations are shown in Fig. 6. It is easy to 
see that in all three cases, the shape of the time dependence of the autocorrelation function 
is the same, and, hence, the slope of the power spectrum does not depend on the cluster core 
radius used in our calculations. 


CONCLUSIONS 


The rotational phase modulation of a pulsar in a globular star cluster may be attributable 
to the fluctuations in the relativistic time delay of the pulsar signal caused by the gravita¬ 
tional held variations due to the motion of the cluster stars. The effect of the low-frequency 
timing noise produced by this process can be studied by accurately determining the behavior 
of the PAT residuals in the time and/or frequency domain. Here, we derived an analytical 
expression for the Shapiro delay of the radio pulse from a pulsar in a globular cluster by 
applying small aberration corrections, of the order of (u/c), where v is the characteristic 
velocity of the stars in the cluster and c is the speed of light. 

Assuming that the interactions between the cluster stars are uncorrelated for large impact 
parameters, the random process produced by the stochastic Shapiro effect manifests itself 
mainly as the noise of the pulsar spin rate. For small impact parameters, we numerically 
simulated the stochastic Shapiro effect by assuming that the mass distribution of the cluster 
stars was specified by the Salpeter function, the velocity distribution of the cluster stars 
was Maxwellian, and the globular cluster density was described either by the model of an 
isothermal sphere or by the King model. 

Our numerical analysis showed that the autocorrelation function of the noise process 
produced by the Shapiro effect attributable to the gravitational held of the moving cluster 
stars could be fitted by a quadratic time polynomial for both the model of an isothermal 
sphere and the King model at large impact parameters. The slope of the power spectrum is 
n ft —1.8 and depends weakly on the model density distribution of the globular star cluster. 
For small impact parameters, the functional time dependence of the autocorrelation function 
for the stochastic Shapiro effect is nearly logarithmic and the slope of the power spectrum 
decreases and is n ft —1.5. 
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Fig. 2: Autocorrelation function of the stochastic Shapiro effect for the model of an 
isothermal sphere with a pulsar at the cluster center for observing times of 1 (solid 
line), 3 (dotted line), and 5 yr (dashed line) 
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Fig. 4: Autocorrelation function of the stochastic Shapiro effect for impact parameters 
smaller than 1 AU in the model of an isothermal sphere with a pulsar at the center 
of a globular cluster for observing times of 1 (solid line), 3 (dotted line), and 5 yr 
(dashed line). 
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Fig. 5: Slope of the power spectrum for the stochastic Shapiro effect for impact pa¬ 
rameters smaller than 1 AU (r = 5 yr). 
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Fig. 6: Autocorrelation function of the stochastic Shapiro effect for the model of an 
isothermal sphere with various globular cluster core radii, 0.52 (1), 0.08 (2), and 0.1 
pc (3), and with a pulsar at the cluster center. The observing time is t = 1.5 yr. 






